
import ij.ImagePlus;
import ij.ImageStack;
import ij.io.FileSaver;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import java.util.Timer;
import static palmer.PALMER_.zoom;

/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

/**
 *
 * @author cloudnine
 */
public class Maliang_runner_ {
    
    //basic 
    public static ImagePlus img_plus = Maliang_initial_.img_plus;
    public static ImageProcessor img_prpcessor = Maliang_initial_.img_prpcessor;
    public static ImageStack img_stack = Maliang_initial_.img_stack;

       
    //malloc CONV_ 卷积参数 
    public static int CONV_KernelLength = (int)Maliang_initial_.OPTION_PSFWinLength + 2;
    public static float [] CONV_Kernel  = new float[CONV_KernelLength];
    public static float [] CONV_Kernel1 = new float[CONV_KernelLength];
    public static float [] CONV_Kernel2 = new float[CONV_KernelLength];
    public static float [] CONV_Kernel3 = new float[CONV_KernelLength];

    //该函数送入GPU处理的参数 HANDLE_
    public static int HANDLE_ConvW = 0, HANDLE_ConvH = 0;

    //location LOC_ 参数
    public static double LOC_threshold;
    public static int LOC_moleculeNum;
    public static int LOC_subregionSize = 7;
    
    //显示 DISPLAY_
    //放大倍数 四舍五入 因为java直接舍掉小数
    public static int DISPLAY_enlarge = (int)(Maliang_initial_.OPTION_PixelSizeRawImage / Maliang_initial_.OPTION_PixelSizeRenderedImage + 0.5) ;
    
    public static void statrRun()
    {    
        
        /*
        -0.0625   -0.0625   -0.0625   -0.0625   -0.0625
        -0.0625    0.1111    0.1111    0.1111   -0.0625
        -0.0625    0.1111    0.1111    0.1111   -0.0625
        -0.0625    0.1111    0.1111    0.1111   -0.0625
        -0.0625   -0.0625   -0.0625   -0.0625   -0.0625
        *///k1*k2 + k3*k4 最终的卷积模板，关于耗时持怀疑态度？有时间一定要测试一下。
        for (int i=0; i<CONV_KernelLength; i++)
        {
                CONV_Kernel[i] = -1f/(float)(CONV_KernelLength*CONV_KernelLength-(CONV_KernelLength-2)*(CONV_KernelLength-2));
                CONV_Kernel1[i] = 1f;
                CONV_Kernel2[i] = 1f/(float)(CONV_KernelLength*CONV_KernelLength-(CONV_KernelLength-2)*(CONV_KernelLength-2))+1f/(float)((CONV_KernelLength-2)*(CONV_KernelLength-2));
                CONV_Kernel3[i] = 1f;	
        }
        CONV_Kernel2[0] = 0; CONV_Kernel2[CONV_KernelLength-1] = 0;
        CONV_Kernel3[0] = 0; CONV_Kernel3[CONV_KernelLength-1] = 0;  
        
        //提取局部最大值的阈值
        LOC_threshold = Maliang_initial_.OPTION_BgThreshold * backNoiseStd(20);
        
        //多帧图片合成一张作为待卷积图像(同时卷积多帧图片,由用户设置一次多少帧)
	HANDLE_ConvW = (Maliang_initial_.IMAGE_imageW + 63) / 64 * 64;
	HANDLE_ConvH = (Maliang_initial_.IMAGE_imageH * Maliang_initial_.OPTION_noFramePerCycle + 15) / 16 * 16;
        
        CUDADLL_.CudaMalloc(HANDLE_ConvW, HANDLE_ConvH, CONV_KernelLength, CONV_Kernel,CONV_Kernel1, CONV_Kernel2, CONV_Kernel3);
    
        // raw PALM image data / de-noising image data
        float[] PALM_rawimage = new float[HANDLE_ConvW * HANDLE_ConvH];
        float[] DeNoise_imageConv = new float[HANDLE_ConvW * HANDLE_ConvH];
        
        //uint8 灰度图像 多帧tif文件保存在一个图像文件中
        if (img_prpcessor instanceof ByteProcessor)
        {
            byte[] imageDataTemp = new byte[Maliang_initial_.IMAGE_imageW * Maliang_initial_.IMAGE_imageH];
            for (int i = 1; i <= Maliang_initial_.OPTION_noFramePerCycle; i++) {
                    imageDataTemp = (byte[]) img_stack.getPixels(i);
                    for (int c = 0; c < Maliang_initial_.IMAGE_imageH; c++) {
                            for (int r = 0; r< Maliang_initial_.IMAGE_imageW; r++) {
                                    PALM_rawimage[(i - 1) * HANDLE_ConvW * Maliang_initial_.IMAGE_imageH + c * HANDLE_ConvW + r] =  (imageDataTemp[c * Maliang_initial_.IMAGE_imageW + r] & 0xFF);    
                            }
                    }
            }
        }
        //saveImage(PALM_rawimage, "子区域提取过程-多帧合成的一张图片", HANDLE_ConvW, HANDLE_ConvH);
        
        // 图像去噪
        CUDADLL_.imageConv(DeNoise_imageConv, PALM_rawimage, HANDLE_ConvW, HANDLE_ConvH, CONV_KernelLength);
        //saveImage(DeNoise_imageConv, "子区域提取过程-去噪结果", HANDLE_ConvW, HANDLE_ConvH);
        
        LOC_moleculeNum = extractSubregion(DeNoise_imageConv, Maliang_initial_.OPTION_noFramePerCycle, LOC_threshold);       
        //saveImage(DeNoise_imageConv, "子区域提取-局部极大值", HANDLE_ConvW, HANDLE_ConvH);
        
        if(LOC_moleculeNum <= 0 ) return;
        
        //single moleucle infromation
        float moleculeInf[] = new float[7 * LOC_moleculeNum]; 
        int[] moleculeLocationPixel = new int[LOC_moleculeNum];
        for (int i = 0; i < LOC_moleculeNum; i++) {
                moleculeLocationPixel[i] = (int) DeNoise_imageConv[i];
        }
        
        //????
	float[] A = { -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, -1,
			0, 0, 0, 0, 0, -1 };
	float[] eye = { 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
			0, 0, 0, 0, 1 };
        
        CUDADLL_.moleculeLoaction(moleculeLocationPixel, LOC_moleculeNum, A,eye, moleculeInf, HANDLE_ConvW, (float)(Maliang_initial_.OPTION_GrayToPhoton * Maliang_initial_.OPTION_QE), LOC_subregionSize);
        
        //moleculeInf 
        // ？   ？    y   x  ？ ？ ？     
        
         // 得到重建真实坐标(基于原始图尺寸)。子区域顶点位置 + 定位位置  ?为什么还减1？是因为定位的时候 不是从0开始 是从1开始的是吧
        for (int i = 0; i < LOC_moleculeNum; i++) {
                moleculeInf[7 * i + 2] += moleculeLocationPixel[i] / HANDLE_ConvW % Maliang_initial_.IMAGE_imageH - (LOC_subregionSize-1)/2 - 1;
                moleculeInf[7 * i + 3] += moleculeLocationPixel[i] % HANDLE_ConvW - (LOC_subregionSize-1)/2 - 1;
        } 
        //saveImage(moleculeInf, "子区域定位-7参数表", 7, LOC_moleculeNum);

        //图像重建 的一些参数
        FloatProcessor fpRebuild = new FloatProcessor(Maliang_initial_.IMAGE_imageW * DISPLAY_enlarge, Maliang_initial_.IMAGE_imageH * DISPLAY_enlarge);
        fpRebuild.setValue(0.0);
        //fpRebuild.setColorModel(hotCM);
        ImagePlus impRebuild = new ImagePlus("重建渲染图", fpRebuild);

        for(int i=0; i < LOC_moleculeNum; i++)
        { 
            int ntemp = 1;
            double gaussValue[][] = new double[DISPLAY_enlarge + 1][DISPLAY_enlarge + 1];
            double gaussValueSum = 0;
            for (int j = -DISPLAY_enlarge / 2; j <= DISPLAY_enlarge / 2; j++) {
                    for (int k = -DISPLAY_enlarge / 2; k <= DISPLAY_enlarge / 2; k++) {
                            gaussValue[k + DISPLAY_enlarge / 2][j + DISPLAY_enlarge / 2] = Math.exp(-1 * (k * k + j * j)/ (2 * ntemp * ntemp));
                            gaussValueSum = gaussValueSum + gaussValue[k + DISPLAY_enlarge / 2][j + DISPLAY_enlarge / 2];
                    }
            }
            
            for (int j = -DISPLAY_enlarge / 2; j <= DISPLAY_enlarge / 2; j++) {
                for (int k = -DISPLAY_enlarge / 2; k <= DISPLAY_enlarge / 2; k++) {
                        fpRebuild.setf((int) (moleculeInf[i * 7 + 3] * DISPLAY_enlarge)+ k,(int) (moleculeInf[i * 7 + 2] * DISPLAY_enlarge)+ j, 
                                (Float.intBitsToFloat(fpRebuild.getPixel((int) (moleculeInf[i * 7 + 3] * DISPLAY_enlarge)+ k, (int) (moleculeInf[i * 7 + 2] * DISPLAY_enlarge)+ j))) + (float) (gaussValue[k+ DISPLAY_enlarge / 2][j + DISPLAY_enlarge / 2] / gaussValueSum));
                }
            }
        }
            impRebuild.resetDisplayRange();
            impRebuild.setDisplayRange(0.0D, fpRebuild.getMax());
            impRebuild.updateAndDraw();
            impRebuild.show();
        
    }
    

    // compute the background deviation 计算图像背景噪声强度标准差
    private static double backNoiseStd(int inFrame)
    {
        //多帧图片合成一张作为待卷积图像(同时卷积多帧图片)
	HANDLE_ConvW = (Maliang_initial_.IMAGE_imageW + 63) / 64 * 64;
	HANDLE_ConvH = (Maliang_initial_.IMAGE_imageH * inFrame + 15) / 16 * 16;
        
        CUDADLL_.CudaMalloc(HANDLE_ConvW, HANDLE_ConvH, CONV_KernelLength, CONV_Kernel,CONV_Kernel1, CONV_Kernel2, CONV_Kernel3);
        
        //原始图像 卷积结果图像
        float[] HANDLE_rawimage  = new float[HANDLE_ConvW * HANDLE_ConvH];
        float[] HANDLE_convimage = new float[HANDLE_ConvW * HANDLE_ConvH];
    
        //uint8 灰度图像 多帧tif文件保存在一个图像文件中
        //note:java中byte是有符号的，java天生没有 unsigned
        //&0xFF 转换为无符号类型 why?
        if (img_prpcessor instanceof ByteProcessor)
        {
            byte[] imageDataTemp = new byte[Maliang_initial_.IMAGE_imageW * Maliang_initial_.IMAGE_imageH];
            for (int i = 1; i <= inFrame; i++) {
                    imageDataTemp = (byte[]) img_stack.getPixels(i);
                    for (int c = 0; c < Maliang_initial_.IMAGE_imageH; c++) {
                            for (int r = 0; r< Maliang_initial_.IMAGE_imageW; r++) {
                                    HANDLE_rawimage[(i - 1) * HANDLE_ConvW * Maliang_initial_.IMAGE_imageH + c * HANDLE_ConvW + r] =  (imageDataTemp[c * Maliang_initial_.IMAGE_imageW + r] & 0xFF);    
                            }
                    }
            }
        }
     
        CUDADLL_.imageConv(HANDLE_convimage, HANDLE_rawimage, HANDLE_ConvW, HANDLE_ConvH, CONV_KernelLength);
      
        //卷积后图像的背景标准差
        double backgroundStd = 0;
        for(int i = CONV_KernelLength; i < Maliang_initial_.IMAGE_imageH - CONV_KernelLength; i++){
            for(int j = CONV_KernelLength; j <  Maliang_initial_.IMAGE_imageW - CONV_KernelLength; j++){
                double maenTemp = 0;
                double stdTemp = 0;
                for(int imgIndex = 0; imgIndex < inFrame; imgIndex++){
                    maenTemp +=  HANDLE_convimage[(imgIndex * HANDLE_ConvW * Maliang_initial_.IMAGE_imageH) + i * HANDLE_ConvW + j];
                }
                //某像素点的均值
                maenTemp = maenTemp / inFrame;
                
                for(int imgIndex = 0; imgIndex < inFrame; imgIndex++){
                    stdTemp +=  Math.pow(HANDLE_convimage[(imgIndex * HANDLE_ConvW * Maliang_initial_.IMAGE_imageH) + i * HANDLE_ConvW + j] - maenTemp, 2);
                }
                stdTemp  = stdTemp / inFrame;
                stdTemp  = Math.sqrt(stdTemp);
                backgroundStd += (float)stdTemp;
        }
        }
        
        backgroundStd = backgroundStd / ((Maliang_initial_.IMAGE_imageW - 2 * CONV_KernelLength) * (Maliang_initial_.IMAGE_imageH - 2 * CONV_KernelLength));
        
        CUDADLL_.CudaFree();

        //saveImage(HANDLE_rawimage, "背景计算-多帧合成的一张图片", HANDLE_ConvW, HANDLE_ConvH);
        //saveImage(HANDLE_convimage, "背景计算-卷积结果", HANDLE_ConvW, HANDLE_ConvH);
        
        return backgroundStd;
    }
    

    
    public static int extractSubregion(float[] image, int inFrame, double threshold)
    {
       int MoleculeCount = 0;

       for(int imgIndex = 0; imgIndex < inFrame; imgIndex++){
           int BasePixel = imgIndex * HANDLE_ConvW * Maliang_initial_.IMAGE_imageH;
           //卷积核半径为2 长度为2*2+1
           //子区域半径尺寸
           int KernelLength = (int)Maliang_initial_.OPTION_PSFWinLength;
           int SubRadiu = (int)Maliang_initial_.OPTION_PSFWinLength - 2; //7*7子区域
           
           for(int i = KernelLength; i < Maliang_initial_.IMAGE_imageH - KernelLength; i++){
               for(int j = KernelLength; j < Maliang_initial_.IMAGE_imageW - KernelLength; j++){

                   if(image[BasePixel + (i * HANDLE_ConvW + j)] >= threshold)
                   {
                       //记录当前点信息
                       float pixelMax = image[BasePixel + (i * HANDLE_ConvW + j)];
                       int x = j;
                       int y = i;

                       //找到该点为7*7邻域中最大值
                       for(int ii = -SubRadiu; ii <= SubRadiu; ii++){
                           for(int jj = -SubRadiu; jj <= SubRadiu; jj++){
                               if(image[BasePixel + (i + ii) * HANDLE_ConvW + (j + jj)] > pixelMax){
                                   pixelMax = image[BasePixel + (i + ii) * HANDLE_ConvW + (j + jj)];
                                   x = (j + jj);
                                   y = (i + ii);
                               }
                           }
                       }

                       //判断最大点是否在去除卷积模板后的区间内
                       if( KernelLength <= x && x < Maliang_initial_.IMAGE_imageW - KernelLength && KernelLength <= y && y < Maliang_initial_.IMAGE_imageH - KernelLength){
                           //位置信息保存在ImageConv头部
                           image[MoleculeCount] = BasePixel + y * HANDLE_ConvW + x;
                           MoleculeCount++;

                           for(int ii = y-SubRadiu; ii <= y+SubRadiu; ii++){
                               for(int jj = x-SubRadiu; jj <= x+SubRadiu; jj++){
                                 image[BasePixel + ii * HANDLE_ConvW + jj] = 0;
                               }
                           }                           
                           image[BasePixel + y * HANDLE_ConvW + x] = 255;
                       }  
                   }
               }
           }
       }
        return MoleculeCount;
    }
    
 
    
    public static void saveImage(float[] Image, String filename, int imageW, int imageH){
        FloatProcessor ip = new FloatProcessor(imageW,imageH,Image,null);
        ImagePlus imp   = new ImagePlus(filename, ip);
        FileSaver file = new FileSaver(imp);
        String saveFileImage = String.format(".\\%s.tif",filename);
        file.saveAsTiff(saveFileImage);
        imp.setDisplayRange(img_plus.getDisplayRangeMin(), img_plus.getDisplayRangeMax());
        imp.show();
    }  
}



